Lecture 9: Extending the Linear Model

Analysing Data

Overview

• Week 8: The Linear Model
• TODAY: The Linear Model with multiple predictors
• [3 weeks for Spring vacation]
• Week 10: Effect sizes
• Week 11: Consolidation

Linear Models (LM): recap

• Trying to make predictions about the world
• Capture the relationship between predictor & outcome
• Linear model (line) is described by an intercept (b₀) and a slope (b₁)

\[Outcome = b_0 + b_1\times Predictor_1 + \varepsilon\]

• Use the model to see:
  • R²: proportion of variance in outcome that model explains
  • t-statistic and associated p-value: is b₁ different from 0?
  • Direction of relationship between predictor & outcome

Today’s Topics

• How good is our model?
  • Our model vs the mean
  • Error in the model
• Throwing more predictors into the mix!
• Comparing hierarchical models
• Comparing predictors in a model
• Bringing it all together: what can our model tell us?

How good is our model?

• A good linear model (LM):
  • Explains a lot about our outcome
  • Captures more than using the simplest model possible
  • Doesn’t contain much error in prediction

Keep is simple, stupid: The mean model

• We want a LM that explains more than the simplest model possible
• The simplest model is the mean…
  • Mean chocolate bars eaten per week is 3.27
    
  • Predict how much chocolate your neighbour eats a week…
• Let’s plot this
  • Chocolate eaten is our outcome variable
  • Not a great model!

• Let’s add a predictor
• Error in model
  • between prediction (using mean) and observed data

• Let’s add a predictor
• Error in model
  • between prediction (using mean) and observed data

We can do better, probably…

• Let’s construct a LM between chocolate liking and eating
  • The line explaining our data with the least error
  • Compare to the line representing our mean model
  • Shows the improvement in prediction from LM (vs mean)

• Let’s construct a LM between chocolate liking and eating
  • The line explaining our data with the least error
  • Compare to the line representing our mean model
  • Shows the improvement in prediction from fitting a LM

• Larger difference = greater improvement!
• We want a LM that explains more than the mean model

Interim summary

• Is our LM is better than the simplest model possible?
  • The mean model is the simplest model
    
  • We want a LM that explains more than the mean model
• But…a LM better than the mean still contains error
• How much error is okay?
  • Does the model explain more than it doesn’t explain?
  • How well does the LM predict the outcome?

Error in the model

• No model will fit the data perfectly
• Fit a line to best capture relationships between variables
• Want the least error possible
  • Compare predicted and observed data points

 

Bringing it all together: the F-statistic

• This is the F-statistic!

\[F =\frac{what\ the model\ can\ explain}{what\ the\ model\ cannot\ explain }=\frac{signal}{noise}\]

• We want signal to be as big as possible

• We want noise to be as small as possible

• A ratio of variance explained relative to varience unexplained

• Ratio > 1 means our model can explain more than it cannot explain

• Associated p-value of how likely we are to find a F-statistic as large as the observed if the null hypothesis is true

Why you gotta complicate things?

 

• Multiple predictors are not much more complicated!
  • We build models to predict what is happening in the world
• Simple explainations for complex relationships?
  
• Multiple predictors = greater explanatory power

How multiple predictors (don’t really) change the LM equation

\[\begin{aligned}Outcome &= Model + Error\\ Y&=b_0 + b_1\times Predictor_1 + \varepsilon\\ &=b_0 + b_1\times Predictor_1 + b_2\times Predictor_2 + \varepsilon\end{aligned}\]

• Y: outcome
  
• b₀: value of outcome when predictors are 0 (the intercept)
  
• b₁: change in outcome associated with a unit change in predictor 1
• b₂: change in outcome associated with a unit change in predictor 2

How multiple predictors (don’t really) change the LM equation

• One predictor LM = regression line
• Two+ predictor LM = regression plane

Puppies, puppies everywhere!

The outcome:

Predictor 1:

Predictor 2:

\[happiness\ =\ b_0 + b_1\times puppies\ +\ b_2\times fluffy + \varepsilon\]

• Y: happiness
  
• b₀: value of happiness when puppies and fluffy are 0 (the intercept)
  
• b₁: change in happiness associated with a unit change in puppies when fluffy is 0
• b₂: change in happiness associated with a unit change in fluffy when puppies is 0

Model fit: F-statistic

## fit the linear model
eg4.mod1.lm <- lm(happiness ~ puppies, data = eg4.data.tib)
## get fit statistics
broom::glance(eg4.mod1.lm)
# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
1     0.744         0.730  1.63      52.4 9.90e-7     2  -37.1  80.2  83.2
# ... with 2 more variables: deviance <dbl>, df.residual <int>

• The F-statistic is 52.38
• The associated p-value is < .001
  • Adding predictors significantly improved model fit
  • Explained significantly more variance than could not explain

Model fit: R² value

eg4.mod1.lm <- lm(happiness ~ puppies, data = eg4.data.tib)

broom::glance(eg4.mod1.lm)
# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
1     0.744         0.730  1.63      52.4 9.90e-7     2  -37.1  80.2  83.2
# ... with 2 more variables: deviance <dbl>, df.residual <int>

• The R² is 0.74
  • 74% of the variance in happiness was explained by puppies and fluffy ratings
• This is just in our observed data!
• Our adjusted R² value was 0.73
  • If we used the same model with the population, we should be able to explain 73% of the variance in happiness

Interim summary, mark II

• The LM can be expanded to include additional predictors
• The model is still described by an intercept (b₀)
• It now includes slopes (bs) for each predictor
• This creates a regression plane instead of a line
• We can assess how good this model is
  • F-ratio and associated p-value
  • R²
• What if we want to compare two linear models?

Comparing linear models

• We can compare models with different numbers of predictors
  • See which model better captures our outcome
  • The models must be ‘hierarchical’
    • 2^nd model has the same predictors as the 1^st model plus extra
    • 3^rd model has the same predictors as the 2^nd model plus extra

Even more puppies!

Predictor 3:

Model 1:

\[\begin{aligned}happiness\ = &\ b_0\ +\\&\ b_1\times puppies\ +\\ &\ b_2\times fluffy + \varepsilon\end{aligned}\]

Model 2:

\[\begin{aligned}happiness\ = &\ b_0\ +\\&\ b_1\times puppies\ +\\ &\ b_2\times fluffy\ +\\ &\ b_3\times dirt + \varepsilon\end{aligned}\]

## Model 1
eg4.mod1.lm <- lm(happiness ~ puppies, data = eg4.data.tib)

## Model 2
eg4.mod2.lm <- lm(happiness ~ puppies + dirt + wet, data = eg4.data.tib)

• First, which is the better model?
• Second, how useful is the best model?

With great (predictive) power comes great responsibility

• Need to assess which is the best model
• Want to compare Model 1 (2 predictors) and Model 2 (3 predictors)
  • Do extra predictors improve the model?
  • Look at R²-change and F-change

R²-change

Model 1:

broom::glance(eg4.mod1.lm)
# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
1     0.744         0.730  1.63      52.4 9.90e-7     2  -37.1  80.2  83.2
# ... with 2 more variables: deviance <dbl>, df.residual <int>

 

Model 2:

broom::glance(eg4.mod2.lm)
# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
1     0.971         0.965 0.583      178. 1.68e-12     4  -15.3  40.7  45.7
# ... with 2 more variables: deviance <dbl>, df.residual <int>

• The R² for Model 2 is 0.97
  • 97% of the variance in happiness is explained by puppies and fluffy).
• We can calculate the change in R² between Model 1 and Model 2
  • 0.97 - 0.74 = 0.23
  • Model 2 explains 23% more variance than Model 1
  • We accounted for 23% more variance in happinessby adding dirt to the model

F-change

eg4.anova <- anova(eg4.mod1.lm,eg4.mod2.lm) %>% broom::tidy()
eg4.anova
# A tibble: 2 x 6
  res.df   rss    df sumsq statistic       p.value
   <dbl> <dbl> <dbl> <dbl>     <dbl>         <dbl>
1     18 47.8     NA  NA        NA   NA           
2     16  5.43     2  42.4      62.4  0.0000000278

• Does Model 2 explain more variance compared to Model 1?
  • Does adding more predictors significantly improve model fit?
• The F-change statistic is 62.39 with a p-value of < .001
• Model 2 was significantly better at explaining variability in happiness than Model 1

Relative contribution of predictor variables

• Now we know Model 2 is the best model…
• Which predictors explain most variance in the outcome variable?
  • Assess the relative contribution of our predictor variables
  • Look at the b-values and associated p-values for each predictor

Number of puppies…

broom::tidy(eg4.mod2.lm, conf.int = TRUE)
# A tibble: 4 x 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)    3.74     0.438       8.55 2.33e- 7   2.81       4.67 
2 puppies        0.995    0.0566     17.6  7.01e-12   0.875      1.11 
3 dirt           0.181    0.0744      2.43 2.75e- 2   0.0228     0.338
4 wet           -0.796    0.0726    -11.0  7.59e- 9  -0.950     -0.642

• puppies significantly predicted happiness
  • b₁ is positive: as puppies increases, happiness increases
  • For every unit increase in puppies, happiness increases by 0.99 units
  • When all other predictors are held constant!

Fluffiness of puppies…

broom::tidy(eg4.mod2.lm, conf.int = TRUE)
# A tibble: 4 x 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)    3.74     0.438       8.55 2.33e- 7   2.81       4.67 
2 puppies        0.995    0.0566     17.6  7.01e-12   0.875      1.11 
3 dirt           0.181    0.0744      2.43 2.75e- 2   0.0228     0.338
4 wet           -0.796    0.0726    -11.0  7.59e- 9  -0.950     -0.642

• fluffy significantly predicted happiness
  • b₂ is fluffy is a positive: as fluffy increases, happiness increases
  • For every unit increase in fluffy, happiness increases by 0.18 units

Dirtiness of puppies…

broom::tidy(eg4.mod2.lm, conf.int = TRUE)
# A tibble: 4 x 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)    3.74     0.438       8.55 2.33e- 7   2.81       4.67 
2 puppies        0.995    0.0566     17.6  7.01e-12   0.875      1.11 
3 dirt           0.181    0.0744      2.43 2.75e- 2   0.0228     0.338
4 wet           -0.796    0.0726    -11.0  7.59e- 9  -0.950     -0.642

• dirt significantly predicted happiness
  • b₃ is dirt is a negative: as dirt increases, happiness decreases
  • For every unit increase in dirt, happiness decreases by -0.8 units

Relative contribution of predictor variables

• We CANNOT interpret the relative contribution from these bs!
• We can’t say that puppies caused happiness to change by a greater amount than dirt
  • puppies and fluffy were measured in different units
  • We can only compare predictors measured in the same unit (e.g. cm or £)
• Instead, need to look at standardized versions of our betas
  • These are ‘measured’ in standard deviations
• Unstandardized b₁ = number of units that happiness increases for every unit increase in puppies
• Standardized b₁ = number of standard deviations (SDs) that happiness increases for every standard deviation increase in puppies

lm.beta::lm.beta(eg4.mod2.lm)

Call:
lm(formula = happiness ~ puppies + dirt + wet, data = eg4.data.tib)

Standardized Coefficients::
(Intercept)     puppies        dirt         wet 
  0.0000000   0.8326795   0.1149789  -0.4675909 

• The standardized beta for puppies (standardized b₁) is 0.83
  • As puppies increase by 1 SD, happiness increases by 0.83 SDs
• The standardized beta for fluffy (standardized b₂) is 0.11
  • As fluffy increases by 1 SD, happiness increases by 0.11 SDs

lm.beta::lm.beta(eg4.mod2.lm)

Call:
lm(formula = happiness ~ puppies + dirt + wet, data = eg4.data.tib)

Standardized Coefficients::
(Intercept)     puppies        dirt         wet 
  0.0000000   0.8326795   0.1149789  -0.4675909 

• The standardized beta for dirt (standardized b₂) is -0.47
  • As dirt increases by 1 SD, happiness decreases by -0.47 SDs

Relative contribution of predictors

lm.beta::lm.beta(eg4.mod2.lm)

Call:
lm(formula = happiness ~ puppies + dirt + wet, data = eg4.data.tib)

Standardized Coefficients::
(Intercept)     puppies        dirt         wet 
  0.0000000   0.8326795   0.1149789  -0.4675909 

• All predictors (bs) are now in the same unit - standard deviations - Now we can directly compare their relative contribution
• 0.83 is the largest standardized beta
• So each SD increase in puppies predicts a larger SD change in happiness than the other predictors

Summary

• Run hierarchical models
• Assess if either model is better than the mean model:
  • F-statistic and associated p-value
• Assess which is better model
  • R²-change, F-change statistic and associated p-value
• Assess which variables significantly predict the outcome
  • t-value and associated p-value
• Assess which variables have largest contribution to the outcome
  • Standardized bs
• The LM captures the relationship between at one or more predictors, x, and an outcome, y
  • LM equation: Outcome = b₀ + b₁ × Predictor 1 + b₂× Predictor 2 + b_n × Predictor n… + Error
• Is our model useful? Assess model fit:
  • F-statistic (or F-change) and associated p-value
  • R² (or R²-change)
• What does our model tell us?
  • Unstandardized bs, t-statistic and associated p-value
  • Standardized bs

Don’t forget to save this page on your computer, otherwise you will lose all notes!